####
# Replication code for "The effect of public safety power shutoffs on climate attitudes and behavioral intentions"
# Matto Mildenberger, Peter D. Howe, Samuel Trachtman, Leah C. Stokes, and Mark Lubell
# Nature Energy
####

rm(list=ls())
setwd("C:/Users/samtr/Dropbox/ENVENT/electricity_outages/replication")

#libaries
library(rio)
library(Matching)
library(MatchIt)
library(gridExtra)
library(tidyverse)
library(estimatr)
library(ggplot2)
library(lmtest)
library(AER)
library(ivpack)
library(plyr)

#survey data
df <- import("outage_survey.csv", header=TRUE)


##### descriptive stats on reason for outage ##### 
outage_reasons <- import("reasons_for_outages.csv")
sum <- nrow(outage_reasons)
sum(outage_reasons$`Fire, wind, or weather`,na.rm=T)/sum
sum(outage_reasons$`PG&E negligence or corruption`,na.rm=T)/sum
sum(outage_reasons$`PG&E safety concerns`,na.rm=T)/sum
sum(outage_reasons$`Government negligence or corruption`,na.rm=T)/sum
sum(outage_reasons$`Climate change`,na.rm=T)/sum


##### figure 2 code #####

bartextsize<-4
logogreen<-"#19a883"
darkergreen<-"#168167"
darkestgreen<-'#0e5141'
almostblack<-'#262626'
midgray<-'#7f7f7f'
lightgray<-'#a6a6a6'
darkgray<-'#595959'
lightestgray<-'#d9d9d9'
marigold<-'#eeac2c'  
darkermarigold<-"#c38d25"
darkestmarigold<-"#6d4f14"
lightmarigold<-"#f4ca78"

inside <- df[which(df$shut_off_planned==1),]
inside$weight <- 1

Q73sub <- inside %>%
  dplyr::select(Q73_1, Q73_2, Q73_3, Q73_4, Q73_5, Q73_6, Q73_7, weight)

Q73samp<-Q73sub%>%
  gather(var,val,-weight)%>%
  mutate(var=as.factor(var))%>%
  group_by(var)%>%
  dplyr::summarise(n=sum(weight))

Q73freq<-Q73sub%>%
  gather(var,val,-weight)%>%
  mutate(var=as.factor(var))%>%
  group_by(var,val)%>%
  dplyr::summarise(freq=sum(weight))%>%
  left_join(Q73samp,by="var")%>%
  mutate(prop=round(freq/n,2))

Q73freq <- na.omit(Q73freq)
Q73freq$var<-ordered(Q73freq$var,levels=c("Q73_1","Q73_2", "Q73_3", "Q73_4", "Q73_5", "Q73_6", "Q73_7"),
                     labels=c("Running out of food", "Running out of water", "Running out of medicine", "Traveling to/from work", "Completing household tasks", "Being able to contact people", "Caring for family"))

Q73freq$val <- ordered(Q73freq$val, levels=c("1","2","3","4"),
                       labels=c("Not at all", "Not very", "Somewhat", "Very"))
Q73freq<-with(Q73freq,Q73freq[order(c(var)),])

Q73freq<-Q73freq%>%
  mutate(per=prop*100)

Q73freq<-ddply(Q73freq, .(var),
               transform, pos = cumsum(prop) - (0.5 * prop)
)
Q73freq$val <- fct_rev(Q73freq$val)

Q73freq<-Q73freq%>%dplyr::filter(per>0)
ggplot(Q73freq, aes(x=var, fill=val, y=prop)) +
  geom_bar(position="fill",stat="identity") +
  geom_text(data=subset(Q73freq,val!="No answer"),position=position_stack(vjust=0.5),aes(label=per),size=bartextsize,color="white")+
  ggtitle("Level of Worry During Outages About..") +
  scale_y_continuous(#breaks=c(0,0.5,1),labels=c("0%","50%","100%"))+
    name="",expand = c(0,0), breaks=c(0,0.5,1), labels=c("0%", "50%", "100%")) +
  scale_x_discrete(name="") +
  scale_fill_manual(values=c("Not at all"=marigold,"Not very" = logogreen,
                             # "No answer"=lightestgray,
                             "Somewhat"=darkergreen, "Very"=darkestgreen))+
  theme(plot.title = element_text(size = 12, face = "bold"), # controls title font + size
        legend.position="top", legend.justification="left", # moves legend to the top and makes it stacked
        legend.title=element_blank(),
        legend.spacing.x = unit(0.1, 'cm'),
        aspect.ratio = 0.5, # of the entire grid (y vs. x axes)
        axis.text.x=element_text(angle=90,size=9, vjust=0), # of the year ticks on the x axis
        axis.text.y=element_text(size=9)) +
  theme( # remove the vertical grid lines
    panel.grid.major.x = element_line(colour="gray"),
    # panel.grid.major.y = element_line(colour = "gray"),
    panel.grid.minor.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"), # add x and y axis lines
    axis.ticks.y=element_blank(),
    axis.ticks.x=element_blank(),
    plot.title=element_text(vjust=-2),
    plot.margin = margin(0, 0.5, 1, 0.25, "cm"))

###########################################################
##### CORE SPECIFICATION-- TREATMENT DETERMINED BY SELF-REPORT #####
############################################################


##### balance without matching-- table 1 ##### 

#exclude socal sample 
nosoc <- subset(df, sample!="socal" & !is.na(shut_off_planned))
treat <- nosoc$shut_off_planned

# covariates
B <- nosoc[,c("pid5", "ideo", "ed", "age", "income", "female", "married", "employed",  "non_english","smoke")] 

Nt = sum(treat==1)
Nc = sum(treat==0)

colnames = c("Var" , "MeanT", "MeanC", "Pval ttest", "Var ratio", "Pval KS", "QQ med diff")
resbal = matrix(NA, nrow=ncol(B)*2 + 2, ncol = 7, dimnames = list(NULL, colnames))
for(i in 1:ncol(B)) {
  
  indx = (!is.na(B[,i]) & treat==1)
  xtr  = B[indx,i]
  indx = (!is.na(B[,i]) & treat==0)
  xco  = B[indx,i]
  
  trt <- c(rep(1, length(xtr)), rep(0,length(xco)))
  var <- c(xtr, xco)
  
  bal<-balanceUV(xtr,xco, ks=TRUE,nboots = 1000, paired=FALSE, match=FALSE)
  
  sigmat<-sqrt(bal$var.Tr)      
  sigmac<-sqrt(bal$var.Co)

  resbal[2*(i-1)+1,1] <- names(B)[i]
  resbal[2*(i-1)+1,2] <- round(bal$mean.Tr,digits=3)
  resbal[2*(i-1)+2,2] <- paste("(",round(sigmat/sqrt(Nt),digits=3),")",sep="")
  resbal[2*(i-1)+1,3] <- round(bal$mean.Co,digits=3)
  resbal[2*(i-1)+2,3] <- paste("(",round(sigmac/sqrt(Nc),digits=3),")",sep="")   
  resbal[2*(i-1)+1,4] <- round(bal$p.value,digits=3)
  resbal[2*(i-1)+1,5] <- round(bal$var.ratio,digits=3)
  resbal[2*(i-1)+1,6] <- round(bal$ks$ks.boot.pvalue,digits=3)  
  resbal[2*(i-1)+1,7] <- round(bal$qqsummary$mediandiff,digits=3)    
}
resbal[ncol(B)*2 + 1, 1] <- "Sample size (T,C)"
resbal[ncol(B)*2 + 1, 2] <- Nt
resbal[ncol(B)*2 + 1, 3] <- Nc
resbal[ncol(B)*2 + 2, 1] <- "Global Balance Test"
resbal[ncol(B)*2 + 2, 2] <- as.numeric(global.test$overall[1])
resbal[ncol(B)*2 + 2, 3] <- as.numeric(global.test$overall[3])
print(resbal) 


##### matching ##### 

#select covariates 
X <- nosoc %>%
  dplyr::select(pid5, married, employed, non_english, kids_in_household, ideo, ed, female, smoke, age, ResponseId, sample, shut_off_planned) %>%
  na.omit

set.seed(12345)
#m.out.nosoc <- matchit(shut_off_planned ~ pid5+ married+ employed+ non_english+ kids_in_household+ ideo+ ed+ female + smoke + age2, data = df.nom,
#                  pop.size=10000, method = "genetic", replace = T, estimand = "ATT") #genetic matching w/ replacement, ATT is default. 

mdtemp <- match.data(m.out.nosoc)
#export(mdtemp, "nosoc.exp.mdtemp.xlsx")
mdtemp <- import("nosoc.exp.mdtemp.xlsx")
m.data <- left_join(subset(mdtemp, select = ResponseId), df, by="ResponseId") #join back to broader df 


##### balance table in matched sample-- table 1 #####

#balance tests by covariate and treatment assessment on survey
treat <- m.data$shut_off_planned

# covariates we would like to have balance on 
B <- m.data[,c("pid5", "ideo", "ed", "age", "income", "female", "married", "employed",  "non_english","smoke")] 

Nt = sum(treat==1)
Nc = sum(treat==0)

colnames = c("Var" , "MeanT", "MeanC", "Pval ttest", "Var ratio", "Pval KS", "QQ med diff")
resbal = matrix(NA, nrow=ncol(B)*2 + 2, ncol = 7, dimnames = list(NULL, colnames))
for(i in 1:ncol(B)) {
  
  indx = (!is.na(B[,i]) & treat==1)
  xtr  = B[indx,i]
  indx = (!is.na(B[,i]) & treat==0)
  xco  = B[indx,i]
  
  trt <- c(rep(1, length(xtr)), rep(0,length(xco)))
  var <- c(xtr, xco)
  
  bal<-balanceUV(xtr,xco, ks=TRUE,nboots = 1000, paired=FALSE, match=FALSE)
  
  sigmat<-sqrt(bal$var.Tr)      
  sigmac<-sqrt(bal$var.Co)

  resbal[2*(i-1)+1,1] <- names(B)[i]
  resbal[2*(i-1)+1,2] <- round(bal$mean.Tr,digits=3)
  resbal[2*(i-1)+2,2] <- paste("(",round(sigmat/sqrt(Nt),digits=3),")",sep="")
  resbal[2*(i-1)+1,3] <- round(bal$mean.Co,digits=3)
  resbal[2*(i-1)+2,3] <- paste("(",round(sigmac/sqrt(Nc),digits=3),")",sep="")   
  resbal[2*(i-1)+1,4] <- round(bal$p.value,digits=3)
  resbal[2*(i-1)+1,5] <- round(bal$var.ratio,digits=3)
  resbal[2*(i-1)+1,6] <- round(bal$ks$ks.boot.pvalue,digits=3)  
  resbal[2*(i-1)+1,7] <- round(bal$qqsummary$mediandiff,digits=3)  
}
resbal[ncol(B)*2 + 1, 1] <- "Sample size (T,C)"
resbal[ncol(B)*2 + 1, 2] <- Nt
resbal[ncol(B)*2 + 1, 3] <- Nc
resbal[ncol(B)*2 + 2, 1] <- "Global Balance Test"
resbal[ncol(B)*2 + 2, 2] <- as.numeric(global.test$overall[1])
resbal[ncol(B)*2 + 2, 3] <- as.numeric(global.test$overall[3])

print(resbal) 
balance.table.matched <- resbal[,c(2:4)]
balance.table.overall <- cbind(balance.table.unmatched, balance.table.matched)


##### policy outcomes-- fig3 #####
outcome.df <- m.data %>% dplyr::select(ces, net2035, worried_warming, shut_off_planned)
outcome.df <- outcome.df %>% dplyr::rename("treat"="shut_off_planned")
cols <- ncol(outcome.df)-1

mean.t = rep(NA, cols)
mean.c = rep(NA, cols)
p.val = rep(NA, cols)
cilow.t <- rep(NA, cols)
cihigh.t <- rep(NA, cols)
cilow.c <- rep(NA,cols)
cihigh.c <- rep(NA, cols)
for (i in 1:cols) {
  mean.t[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[2]
  cilow.t[i] <-  t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[2]
  mean.c[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[1]
  cilow.c[i] <-  t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[1]
  cihigh.c[i] <- t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[2]
  p.val[i] <- t.test(outcome.df[,i]~outcome.df$treat)$p.value
}
means <- c(mean.t, mean.c)
cilow <- c(cilow.t,cilow.c)
cihigh <- c(cihigh.t, cihigh.c)
treat <- c(rep("Treat", cols), rep("Control", cols))
plotdata <- data.frame(rep(names(outcome.df)[1:ncol(outcome.df)-1],2) , means, cilow, cihigh, treat)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='ces', 1,
                              ifelse(outcome=='net2035', 2,
                                     ifelse(outcome=='worried_warming', 3, NA))))

plotdata$order <- (cols+1)-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$zone <- ifelse(plotdata$treat=="Treat", "Outage-exposed respondents", "Matched non-outage respondents")
plotdata$x <- c(.9,1.1,1.9,2.1,2.9,3.1)
outcomes2 <- c("80% clean energy by 2030", "Net-0 carbon by 2035", "Worried about global warming")

policy <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=zone))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=zone))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "")+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))+
  scale_y_continuous(limits=c(1,4), labels=c("Very opposed", "Somewhat opposed", "Somewhat supportive","Very supportive"))+
  scale_color_manual(values=c('grey','black'))+
  theme(legend.title=element_blank(), legend.position="bottom")+
  theme(text = element_text(size=14)) #+
policy
ggsave("figs/policy_experience.png", policy, width=8, height=5)

##### behavioral intentions- fig4 and appendix #####

outcome.df <- m.data %>% dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping, shut_off_planned) 
cols <- ncol(outcome.df)-1
outcome.df <- outcome.df %>% dplyr::rename("treat"="shut_off_planned")

mean.t = rep(NA, cols)
mean.c = rep(NA, cols)
p.val = rep(NA, cols)
cilow.t <- rep(NA, cols)
cihigh.t <- rep(NA, cols)
cilow.c <- rep(NA,cols)
cihigh.c <- rep(NA, cols)
for (i in 1:cols) {
  mean.t[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[2]
  cilow.t[i] <-  t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[2]
  mean.c[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[1]
  cilow.c[i] <-  t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[1]
  cihigh.c[i] <- t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[2]
  p.val[i] <- t.test(outcome.df[,i]~outcome.df$treat)$p.value
}
means <- c(mean.t, mean.c)
cilow <- c(cilow.t,cilow.c)
cihigh <- c(cihigh.t, cihigh.c)
treat <- c(rep("Treat", cols), rep("Control", cols))
plotdata <- data.frame(rep(names(outcome.df)[1:ncol(outcome.df)-1],2) , means, cilow, cihigh, treat)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- c(.9,1.1,1.9,2.1,2.9,3.1,3.9,4.1,4.9,5.1, 5.9,6.1,6.9,7.1)
plotdata$zone <- ifelse(plotdata$treat=="Treat", "Outage-exposed respondents", "Matched non-outage respondents")
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")

adaptation.plot <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=zone))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=zone))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "Proportion intending to take adaptation behavior")+
  ylim(0,.7)+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  scale_color_manual(values=c('grey','black'))+
  theme(legend.title=element_blank())+
  ggplot2::annotate("text", x = 7.2, y = .18, label = "**")+
  ggplot2::annotate("text", x = 3.2, y = .18, label = "**")+
  ggplot2::annotate("text", x = 1.2, y = .22, label = "*")+
  theme(text = element_text(size=14)) +
  theme(legend.title=element_blank(), legend.position="bottom")
adaptation.plot
ggsave("figs/adaptation_experience.png", adaptation.plot, width=7, height=5)

#p-vals
options(scipen=100)
t.test(m.data$backup ~m.data$shut_off_planned)$p.value
t.test(m.data$landscaping ~m.data$shut_off_planned)$p.value
t.test(m.data$ev ~m.data$shut_off_planned)$p.value

lm_robust(backup ~ shut_off_planned, data = m.data)
lm_robust(landscaping ~ shut_off_planned, data = m.data)
lm_robust(ev ~ shut_off_planned, data = m.data)


### adaptation by global warming beliefs-- fig5 
outcome.df.believe <- m.data %>% 
  filter(warming_anthro==1) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping,shut_off_planned)

outcome.df.deny <- m.data %>% 
  filter(warming_anthro==0) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping,shut_off_planned)

cols <- ncol(outcome.df.believe)-1

coef.believe = rep(NA, cols)
se.believe = rep(NA, cols)
for (i in 1:cols) {
  reg <- lm_robust(outcome.df.believe[,i] ~ outcome.df.believe$shut_off_planned)
  coef.believe[i] <- coeftest(reg)[2,1]
  se.believe[i] <- coeftest(reg)[2,2]
}

coef.deny = rep(NA, cols)
se.deny = rep(NA, cols)
for (i in 1:cols) {
  reg <- lm_robust(outcome.df.deny[,i] ~ outcome.df.deny$shut_off_planned)
  coef.deny[i] <- coeftest(reg)[2,1]
  se.deny[i] <- coeftest(reg)[2,2]
}

plotdata <- data.frame(c(coef.believe, coef.deny), c(se.believe, se.deny), rep(names(outcome.df.believe)[1:7], 2), c(rep("Accept",7), rep("Deny", 7)))
names(plotdata) <- c("coef","se","outcome", "group")
plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
perms <- 2
x <- sort(c((1:cols-.1),(1:cols+.1)))
plotdata$x <- x

het_effects_gw <- ggplot(plotdata, aes(x = x, y = coef))+
  geom_point(size=3, aes(shape = group, color = group))+
  geom_linerange(aes(ymin=coef-1.96*se, ymax=coef+1.96*se, color=group), linetype = 1, size = 1)+
  geom_linerange(aes(ymin=coef-1.645*se, ymax=coef+1.645*se, color=group), linetype = 1, size = 1.5)+  coord_flip()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme_classic()+
  labs(x = "", y = "Estimated effect of outage on adaptation intentions")+
  ylim(-.5,.5)+
  scale_x_continuous(breaks = 1:cols, labels=outcomes2)+ #, labels = outcomes2
  theme(text = element_text(size=14))+
  scale_color_manual(values=c('grey','brown'))+
  scale_shape_manual(values=c("triangle","square"))+
  theme(legend.position = "bottom", legend.title=element_blank())
het_effects_gw
ggsave("figs/het_effects_gw.png", het_effects_gw, width=7, height=5)
table(m.data$warming_anthro)

### adaptation-- robustness to covariate adjustment, SI2 
coef = rep(NA, cols)
se = rep(NA, cols)
for (i in 1:cols) {
  reg <- lm_robust(outcome.df[,i] ~ outcome.df$treat + m.data$female + m.data$ideo + m.data$ed + m.data$age + m.data$income + m.data$pid5 + m.data$kids_in_household + 
                     m.data$employed + m.data$non_english)
  coef[i] <- coeftest(reg)[2,1]
  se[i] <- coeftest(reg)[2,2]
}
coef

plotdata <- data.frame(coef, se, names(outcome.df)[1:7])
names(plotdata) <- c("coef","se","outcome")
plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")

cov_adjusted_adpt <- ggplot(plotdata, aes(x = order, y = coef))+
  geom_point(size=3)+
  geom_linerange(aes(ymin=coef-1.96*se, ymax=coef+1.96*se), linetype = 1, size = 1)+
  geom_linerange(aes(ymin=coef-1.645*se, ymax=coef+1.645*se), linetype = 1, size = 1.5)+  coord_flip()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme_classic()+
  labs(x = "", y = "Estimated effect of outage on adaptation intentions")+
  ylim(-.5,.5)+
  scale_x_continuous(breaks = 1:cols, labels=outcomes2)+ #, labels = outcomes2
  theme(text = element_text(size=14))+
  scale_color_manual(values=c('black','grey'))+
  theme(legend.position = "bottom", legend.title=element_blank())
cov_adjusted_adpt
ggsave("figs/cov_adjusted_adpt_experience.png", cov_adjusted_adpt, width=7, height=5)

#pvals
lm_robust(m.data$backup ~ shut_off_planned + female + ideo + ed + age + income + pid5 + kids_in_household + employed + non_english, data = m.data)$p.value
lm_robust(m.data$landscaping ~ shut_off_planned + female + ideo + ed + age + income + pid5 + kids_in_household + employed + non_english, data = m.data)$p.value
lm_robust(m.data$ev ~ shut_off_planned + female + ideo + ed + age + income + pid5 + kids_in_household + employed + non_english, data = m.data)$p.value


### adaptation by income-- SI3
m.data$over100 <- m.data$income >= 3
table(m.data$over100)

outcome.df.over <- m.data %>% 
  filter(over100==1) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping,shut_off_planned)

outcome.df.under <- m.data %>% 
  filter(over100==0) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping,shut_off_planned)

cols <- ncol(outcome.df.over)-1

coef.over = rep(NA, cols)
se.over = rep(NA, cols)
for (i in 1:cols) {
  reg <- lm_robust(outcome.df.over[,i] ~ outcome.df.over$shut_off_planned)
  coef.over[i] <- coeftest(reg)[2,1]
  se.over[i] <- coeftest(reg)[2,2]
}

coef.under = rep(NA, cols)
se.under = rep(NA, cols)
for (i in 1:cols) {
  reg <- lm_robust(outcome.df.under[,i] ~ outcome.df.under$shut_off_planned)
  coef.under[i] <- coeftest(reg)[2,1]
  se.under[i] <- coeftest(reg)[2,2]
}

plotdata <- data.frame(c(coef.over, coef.under), c(se.over, se.under), rep(names(outcome.df.over)[1:7], 2), c(rep("Income over $100,000",7), rep("Income under $100,000", 7)))
names(plotdata) <- c("coef","se","outcome", "group")
plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- x
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")


het_effects.income <- ggplot(plotdata, aes(x = x, y = coef))+
  geom_point(size=3, aes(shape = group, color = group))+
  geom_linerange(aes(ymin=coef-1.96*se, ymax=coef+1.96*se, color=group), linetype = 1, size = 1)+
  geom_linerange(aes(ymin=coef-1.645*se, ymax=coef+1.645*se, color=group), linetype = 1, size = 1.5)+  coord_flip()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme_classic()+
  labs(x = "", y = "Estimated effect of outage on adaptation intentions")+
  ylim(-.5,.5)+
  scale_x_continuous(breaks = 1:cols, labels=outcomes2)+ #, labels = outcomes2
  theme(text = element_text(size=14))+
  scale_color_manual(values=c('grey','black'))+
  scale_shape_manual(values=c("triangle","square"))+
  theme(legend.position = "bottom", legend.title=element_blank())
het_effects.income
ggsave("figs/het_effects_income_experience.png", het_effects.income, width=7, height=5)
table(m.data$over100)

### adaptation by PID-- SI3
m.data$pid <- ifelse(m.data$pid5 < 3, "Democrat",
                     ifelse(m.data$pid5 > 3, "Republican", "Independent"))

outcome.df.dem <- m.data %>% 
  filter(pid=="Democrat") %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping,shut_off_planned) #dropped home materials since so rare. 
outcome.df.rep <- m.data %>% 
  filter(pid=="Republican") %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping,shut_off_planned) #dropped home materials since so rare. 
outcome.df.ind <- m.data %>% 
  filter(pid=="Independent") %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping,shut_off_planned) #dropped home materials since so rare. 


cols <- ncol(outcome.df)-1
coef.dem = rep(NA, cols)
se.dem = rep(NA, cols)
coef.rep = rep(NA, cols)
se.rep= rep(NA, cols)
coef.ind= rep(NA, cols)
se.ind= rep(NA, cols)
for (i in 1:cols) {
  reg.dem <- lm_robust(outcome.df.dem[,i] ~ outcome.df.dem$shut_off_planned)
  coef.dem[i] <- coeftest(reg.dem)[2,1]
  se.dem[i] <- coeftest(reg.dem)[2,2]
  #
  reg.rep <- lm_robust(outcome.df.rep[,i] ~ outcome.df.rep$shut_off_planned)
  coef.rep[i] <- coeftest(reg.rep)[2,1]
  se.rep[i] <- coeftest(reg.rep)[2,2]
  #
  reg.ind <- lm_robust(outcome.df.ind[,i] ~ outcome.df.ind$shut_off_planned)
  coef.ind[i] <- coeftest(reg.ind)[2,1]
  se.ind[i] <- coeftest(reg.ind)[2,2]
}


plotdata <- data.frame(c(coef.dem, coef.rep, coef.ind), c(se.dem, se.rep, se.ind), names(outcome.df)[1:7], c(rep("Democrat", cols), rep("Republican", cols), rep("Independent", cols)))
names(plotdata) <- c("coef","se","outcome", "pid")
plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- c(.9,1, 1.1,1.9,2, 2.1,2.9,3, 3.1,3.9,4, 4.1,4.9,5, 5.1, 5.9,6,6.1,6.9,7, 7.1)
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")

het_effects_pid <- ggplot(plotdata, aes(x = x, y = coef))+
  geom_point(size=3, aes(color=pid, shape = pid))+
  geom_linerange(aes(ymin=coef-1.96*se, ymax=coef+1.96*se, color=pid), linetype = 1, size = 1)+
  geom_linerange(aes(ymin=coef-1.645*se, ymax=coef+1.645*se, color=pid), linetype = 1, size = 1.5)+  coord_flip()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme_classic()+
  labs(x = "", y = "Estimated effect of outage on adaptation intentions")+
  ylim(-.6,.6)+
  scale_x_continuous(breaks = 1:cols, labels=outcomes2)+ #, labels = outcomes2
  theme(text = element_text(size=14))+
  scale_color_manual(values=c('blue','grey', 'red'))+
  theme(legend.position = "bottom", legend.title=element_blank())
het_effects_pid
ggsave("figs/het_effects_pid_experience.png", het_effects_pid, width=7, height=5)
table(m.data$pid)




### adaptation by distance to cutpoint-- SI3
table(m.data$days_outage)
outcome.df.c2 <- m.data %>% 
  filter(shut_off_planned==0 & dist < -1048.94) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)
outcome.df.c1 <- m.data %>% 
  filter(shut_off_planned==0 & dist >= -1048.94) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)
outcome.df.t <- m.data %>% 
  filter(shut_off_planned==1) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)

mean.c1 = rep(NA, cols)
mean.c2 = rep(NA, cols)
mean.t = rep(NA, cols)
cilow.c1 <- rep(NA, cols)
cihigh.c1 <- rep(NA, cols)
cilow.c2 <- rep(NA, cols)
cihigh.c2 <- rep(NA, cols)
cilow.t <- rep(NA,cols)
cihigh.t <- rep(NA, cols)
for (i in 1:cols) {
  mean.c1[i] <- t.test(outcome.df.c1[i])$estimate[1]
  cilow.c1[i] <-  t.test(outcome.df.c1[i])$conf.int[1]
  cihigh.c1[i] <- t.test(outcome.df.c1[i])$conf.int[2]
  mean.c2[i] <- t.test(outcome.df.c2[i])$estimate[1]
  cilow.c2[i] <-  t.test(outcome.df.c2[i])$conf.int[1]
  cihigh.c2[i] <- t.test(outcome.df.c2[i])$conf.int[2]
  mean.t[i] <- t.test(outcome.df.t[i])$estimate[1]
  cilow.t[i] <-  t.test(outcome.df.t[i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df.t[i])$conf.int[2]
}
means <- c(mean.c1, mean.c2, mean.t)
cilow <- c(cilow.c1,cilow.c2, cilow.t)
cihigh <- c(cihigh.c1, cihigh.c2, cihigh.t)
shut_off_planned <- c(rep("Far from outages", cols), rep("Close to outages", cols), rep("Outage-exposed", cols))
plotdata <- data.frame(rep(names(outcome.df.t),3) , means, cilow, cihigh, shut_off_planned)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- c(.9,1, 1.1,1.9,2, 2.1,2.9,3, 3.1,3.9,4, 4.1,4.9,5, 5.1, 5.9,6,6.1,6.9,7, 7.1)
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")

het_adapt_dist <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=shut_off_planned, shape=shut_off_planned))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=shut_off_planned))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "Proportion intending to take adaptation behavior")+
  ylim(-.1,.7)+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  scale_color_manual(values=c('grey','black', 'forest green'))+
  theme(text = element_text(size=14)) +
  theme(legend.position = "bottom", legend.title=element_blank())
het_adapt_dist
ggsave("figs/het_adapt_dist_experience.png", het_adapt_dist, width=7, height=5)


### adaptation by friends exposed-- SI3
outcome.df.c2 <- m.data %>% 
  filter(shut_off_planned==0 & close_friends ==0 ) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)
outcome.df.c1 <- m.data %>% 
  filter(shut_off_planned==0 & close_friends == 1) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)
outcome.df.t <- m.data %>% 
  filter(shut_off_planned==1) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)

mean.c1 = rep(NA, cols)
mean.c2 = rep(NA, cols)
mean.t = rep(NA, cols)
cilow.c1 <- rep(NA, cols)
cihigh.c1 <- rep(NA, cols)
cilow.c2 <- rep(NA, cols)
cihigh.c2 <- rep(NA, cols)
cilow.t <- rep(NA,cols)
cihigh.t <- rep(NA, cols)
for (i in 1:cols) {
  mean.c1[i] <- t.test(outcome.df.c1[i])$estimate[1]
  cilow.c1[i] <-  t.test(outcome.df.c1[i])$conf.int[1]
  cihigh.c1[i] <- t.test(outcome.df.c1[i])$conf.int[2]
  mean.c2[i] <- t.test(outcome.df.c2[i])$estimate[1]
  cilow.c2[i] <-  t.test(outcome.df.c2[i])$conf.int[1]
  cihigh.c2[i] <- t.test(outcome.df.c2[i])$conf.int[2]
  mean.t[i] <- t.test(outcome.df.t[i])$estimate[1]
  cilow.t[i] <-  t.test(outcome.df.t[i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df.t[i])$conf.int[2]
}
means <- c(mean.c1, mean.c2, mean.t)
cilow <- c(cilow.c1,cilow.c2, cilow.t)
cihigh <- c(cihigh.c1, cihigh.c2, cihigh.t)
shut_off_planned <- c(rep("No friends / family outages", cols), rep("Friends / family had outages", cols), rep("Outage-exposed", cols))
plotdata <- data.frame(rep(names(outcome.df.t),3) , means, cilow, cihigh, shut_off_planned)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- c(.9,1, 1.1,1.9,2, 2.1,2.9,3, 3.1,3.9,4, 4.1,4.9,5, 5.1, 5.9,6,6.1,6.9,7, 7.1)
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")

het_adapt_friends <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=shut_off_planned, shape=shut_off_planned))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=shut_off_planned))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "Proportion intending to take adaptation behavior")+
  ylim(-.1,.7)+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  scale_color_manual(values=c('grey','black', 'forest green'))+
  theme(text = element_text(size=14)) +
  theme(legend.position = "bottom", legend.title=element_blank())
het_adapt_friends
ggsave("figs/het_adapt_friends_experience.png", het_adapt_friends, width=7, height=5)


### adaptation by time exposed-- SI3
table(m.data$days_outage)
outcome.df.c <- m.data %>% 
  filter(shut_off_planned==0) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)
outcome.df.t1 <- m.data %>% 
  filter(shut_off_planned==1 & days_outage > 0 & days_outage <= 3) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)
outcome.df.t2 <- m.data %>% 
  filter(shut_off_planned==1  & days_outage > 3) %>%
  dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping)

mean.t1 = rep(NA, cols)
mean.t2 = rep(NA, cols)
mean.c = rep(NA, cols)
cilow.t1 <- rep(NA, cols)
cihigh.t1 <- rep(NA, cols)
cilow.t2 <- rep(NA, cols)
cihigh.t2 <- rep(NA, cols)
cilow.c <- rep(NA,cols)
cihigh.c <- rep(NA, cols)
for (i in 1:cols) {
  mean.t1[i] <- t.test(outcome.df.t1[i])$estimate[1]
  cilow.t1[i] <-  t.test(outcome.df.t1[i])$conf.int[1]
  cihigh.t1[i] <- t.test(outcome.df.t1[i])$conf.int[2]
  mean.t2[i] <- t.test(outcome.df.t2[i])$estimate[1]
  cilow.t2[i] <-  t.test(outcome.df.t2[i])$conf.int[1]
  cihigh.t2[i] <- t.test(outcome.df.t2[i])$conf.int[2]
  mean.c[i] <- t.test(outcome.df.c[i])$estimate[1]
  cilow.c[i] <-  t.test(outcome.df.c[i])$conf.int[1]
  cihigh.c[i] <- t.test(outcome.df.c[i])$conf.int[2]
}
means <- c(mean.t1, mean.t2, mean.c)
cilow <- c(cilow.t1,cilow.t2, cilow.c)
cihigh <- c(cihigh.t1, cihigh.t2, cihigh.c)
shut_off_planned <- c(rep("Less than 3 days", cols), rep("More than 3 days", cols), rep("Matched non-outage", cols))
plotdata <- data.frame(rep(names(outcome.df.c),3) , means, cilow, cihigh, shut_off_planned)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- c(.9,1, 1.1,1.9,2, 2.1,2.9,3, 3.1,3.9,4, 4.1,4.9,5, 5.1, 5.9,6,6.1,6.9,7, 7.1)
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")

het_adapt_time <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=shut_off_planned, shape=shut_off_planned))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=shut_off_planned))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "Proportion intending to take adaptation behavior")+
  ylim(0,.7)+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  scale_color_manual(values=c('grey','black', 'forest green'))+
  theme(text = element_text(size=14)) +
  theme(legend.position = "bottom", legend.title=element_blank())
het_adapt_time
ggsave("figs/het_adapt_time_experience.png", het_adapt_time, width=7, height=5)



##### utility trust-- fig6 ##### 

m.data$utility_trust_binary2 <- 1 - m.data$utility_trust_binary
outcome.df <- m.data %>% dplyr::select(utility_trust_binary2, utility_responsible_binary, pge_liable_binary, pge_role_binary, shut_off_planned) #utility responsibility not as good b/c makes less sense for SoCal 
outcome.df <- outcome.df %>% dplyr::rename("treat"="shut_off_planned")
cols <- ncol(outcome.df)-1

mean.t = rep(NA, cols)
mean.c = rep(NA, cols)
p.val = rep(NA, cols)
cilow.t <- rep(NA, cols)
cihigh.t <- rep(NA, cols)
cilow.c <- rep(NA,cols)
cihigh.c <- rep(NA, cols)
for (i in 1:cols) {
  mean.t[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[2]
  cilow.t[i] <-  t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[2]
  mean.c[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[1]
  cilow.c[i] <-  t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[1]
  cihigh.c[i] <- t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[2]
  p.val[i] <- t.test(outcome.df[,i]~outcome.df$treat)$p.value
}
means <- c(mean.t, mean.c)
cilow <- c(cilow.t,cilow.c)
cihigh <- c(cihigh.t, cihigh.c)
treat <- c(rep("Treat", cols), rep("Control", cols))
plotdata <- data.frame(rep(names(outcome.df)[1:ncol(outcome.df)-1],2) , means, cilow, cihigh, treat)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='utility_trust_binary2', 1,
                              ifelse(outcome=='utility_responsible_binary', 2,
                                     ifelse(outcome=='pge_liable_binary', 3, 
                                            ifelse(outcome=='pge_role_binary', 4, NA)))))

plotdata$order <- (cols+1)-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- c(.9,1.1,1.9,2.1,2.9,3.1, 3.9, 4.1)
plotdata$zone <- ifelse(plotdata$treat=="Treat", "Outage-exposed respondents", "Matched non-outage respondents")
outcomes2 <- c( "PG&E major restructuring", "PG&E liable for wildfires (strongly agree)" , "Utility completely responsible for outages",  "Completely distrust their utility")

utilities <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=zone))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=zone))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "Proportion taking position")+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))+
  scale_color_manual(values=c('grey','black'))+
  theme(legend.position = "bottom", legend.title=element_blank()) +
  ggplot2::annotate("text", x = 4.2, y = .4, label = "**")+
  ggplot2::annotate("text", x = 3.2, y = .71, label = "**")+
  theme(text = element_text(size=14))
utilities
ggsave("figs/utilities_experience.png", utilities, width=7, height=4)
t.test(m.data$utility_trust_binary2 ~ m.data$shut_off_planned)
t.test(m.data$utility_responsible_binary ~ m.data$shut_off_planned)


### utility trust with covariate adjustment-- SI2

coef = rep(NA, cols)
se = rep(NA, cols)
for (i in 1:cols) {
  reg <- lm_robust(outcome.df[,i] ~ outcome.df$treat + m.data$female + m.data$ideo + m.data$ed + m.data$age + m.data$income + m.data$pid5 + m.data$kids_in_household + 
                     m.data$employed + m.data$non_english)
  coef[i] <- coeftest(reg)[2,1]
  se[i] <- coeftest(reg)[2,2]
}


plotdata <- data.frame(coef, se, names(outcome.df)[1:4])
names(plotdata) <- c("coef","se","outcome")
plotdata$order <- with(plotdata,
                       ifelse(outcome=='utility_trust_binary2', 1,
                              ifelse(outcome=='utility_responsible_binary', 2,
                                     ifelse(outcome=='pge_liable_binary', 3, 
                                            ifelse(outcome=='pge_role_binary', 4, NA)))))
plotdata$order <- (cols+1)-plotdata$order
plotdata <- arrange(plotdata, order)

cov_adjusted_utility <- ggplot(plotdata, aes(x = order, y = coef))+
  geom_point(size=3)+
  geom_linerange(aes(ymin=coef-1.96*se, ymax=coef+1.96*se), linetype = 1, size = 1)+
  geom_linerange(aes(ymin=coef-1.645*se, ymax=coef+1.645*se), linetype = 1, size = 1.5)+  coord_flip()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme_classic()+
  labs(x = "", y = "Estimated effect of outage on PG&E attitudes")+
  ylim(-.5,.5)+
  scale_x_continuous(breaks = 1:cols, labels=outcomes2)+ #, labels = outcomes2
  theme(text = element_text(size=14))+
  scale_color_manual(values=c('black','grey'))+
  theme(legend.position = "bottom", legend.title=element_blank())
cov_adjusted_utility
ggsave("figs/cov_adjusted_utility_experience.png", cov_adjusted_utility, width=7, height=5)

#pvals
lm_robust(m.data$utility_trust_binary ~ shut_off_planned + female + ideo + ed + age + income + pid5 + kids_in_household + employed + non_english, data = m.data)$p.value


##### politician approval-- fig7 #####

outcome.df <- m.data %>% dplyr::select(trump_approve, newsom_approve, local_govt_approve, shut_off_planned)
outcome.df <- outcome.df %>% dplyr::rename("treat"="shut_off_planned")
cols <- ncol(outcome.df)-1

mean.t = rep(NA, cols)
mean.c = rep(NA, cols)
p.val = rep(NA, cols)
cilow.t <- rep(NA, cols)
cihigh.t <- rep(NA, cols)
cilow.c <- rep(NA,cols)
cihigh.c <- rep(NA, cols)
for (i in 1:cols) {
  mean.t[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[2]
  cilow.t[i] <-  t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[2]
  mean.c[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[1]
  cilow.c[i] <-  t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[1]
  cihigh.c[i] <- t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[2]
  p.val[i] <- t.test(outcome.df[,i]~outcome.df$treat)$p.value
}
means <- c(mean.t, mean.c)
cilow <- c(cilow.t,cilow.c)
cihigh <- c(cihigh.t, cihigh.c)
treat <- c(rep("Treat", cols), rep("Control", cols))
plotdata <- data.frame(rep(names(outcome.df)[1:ncol(outcome.df)-1],2) , means, cilow, cihigh, treat)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='trump_approve', 1,
                              ifelse(outcome=='newsom_approve', 2,
                                     ifelse(outcome=='local_govt_approve', 3, NA))))

plotdata$order <- (cols+1)-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$zone <- ifelse(plotdata$treat=="Treat", "Outage-exposed respondents", "Matched non-outage respondents")
plotdata$x <- c(.9,1.1,1.9,2.1,2.9,3.1)
outcomes2 <- c( "Local govt approval",  "Newsom approval","Trump approval" )

politics <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=zone))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=zone))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "Approval rating (0-100)")+
  ylim(15,65)+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))+
  scale_color_manual(values=c('grey','black'))+
  ggtitle("(a)")+
  theme(legend.title=element_blank(), legend.position="bottom")+
  theme(text = element_text(size=14)) #+
politics

#split by pid
m.data$pid <- ifelse(m.data$pid5 < 3, "Democrat",
                     ifelse(m.data$pid5 > 3, "Republican", "Independent"))

outcome.df.dem <- m.data %>% 
  filter(pid=="Democrat") %>%
  dplyr::select(trump_approve, newsom_approve, local_govt_approve, shut_off_planned)  
outcome.df.rep <- m.data %>% 
  filter(pid=="Republican") %>%
  dplyr::select(trump_approve, newsom_approve, local_govt_approve, shut_off_planned)  
outcome.df.ind <- m.data %>% 
  filter(pid=="Independent") %>%
  dplyr::select(trump_approve, newsom_approve, local_govt_approve, shut_off_planned)  


cols <- ncol(outcome.df.dem)-1
coef.dem = rep(NA, cols)
se.dem = rep(NA, cols)
coef.rep = rep(NA, cols)
se.rep= rep(NA, cols)
coef.ind= rep(NA, cols)
se.ind= rep(NA, cols)
for (i in 1:cols) {
  reg.dem <- lm_robust(outcome.df.dem[,i] ~ outcome.df.dem$shut_off_planned)
  coef.dem[i] <- coeftest(reg.dem)[2,1]
  se.dem[i] <- coeftest(reg.dem)[2,2]
  #
  reg.rep <- lm_robust(outcome.df.rep[,i] ~ outcome.df.rep$shut_off_planned)
  coef.rep[i] <- coeftest(reg.rep)[2,1]
  se.rep[i] <- coeftest(reg.rep)[2,2]
  #
  reg.ind <- lm_robust(outcome.df.ind[,i] ~ outcome.df.ind$shut_off_planned)
  coef.ind[i] <- coeftest(reg.ind)[2,1]
  se.ind[i] <- coeftest(reg.ind)[2,2]
}


plotdata <- data.frame(c(coef.dem, coef.rep, coef.ind), c(se.dem, se.rep, se.ind), rep(names(outcome.df.dem)[1:3], 3), c(rep("Democrat", cols), rep("Republican", cols), rep("Independent", cols)))
names(plotdata) <- c("coef","se","outcome", "pid")
plotdata$order <- with(plotdata,
                       ifelse(outcome=='trump_approve', 1,
                              ifelse(outcome=='newsom_approve', 2,
                                     ifelse(outcome=='local_govt_approve', 3, NA))))

plotdata$order <- (cols+1)-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- c(.9,1, 1.1,1.9,2, 2.1,2.9,3, 3.1)
outcomes2 <- c( "Local govt approval",  "Newsom approval","Trump approval" )


het_effects_pid <- ggplot(plotdata, aes(x = x, y = coef))+
  geom_point(size=3, aes(color=pid, shape = pid))+
  geom_linerange(aes(ymin=coef-1.96*se, ymax=coef+1.96*se, color=pid), linetype = 1, size = 1)+
  geom_linerange(aes(ymin=coef-1.645*se, ymax=coef+1.645*se, color=pid), linetype = 1, size = 1.5)+  coord_flip()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme_classic()+
  labs(x = "", y = "Estimated effect on approval rating")+
  ylim(-30,30)+
  scale_x_continuous(breaks = 1:cols, labels=outcomes2)+ 
  theme(text = element_text(size=14))+
  scale_color_manual(values=c('blue','grey', 'red'))+
  ggtitle("(b)")+
  theme(legend.position = "bottom", legend.title=element_blank())
het_effects_pid

politics_experience <- grid.arrange(politics, het_effects_pid, ncol=2)
politics_experience
ggsave("figs/politics_experience.png", politics_experience, width=11, height=5)



##### differences in prior adaptation #####

outcome.df <- m.data %>% dplyr::select(already_battery, already_solar, already_backup, already_landscaping, already_ev, shut_off_planned)  
cols <- ncol(outcome.df)-1

mean.t = rep(NA, cols)
mean.c = rep(NA, cols)
p.val = rep(NA, cols)
cilow.t <- rep(NA, cols)
cihigh.t <- rep(NA, cols)
cilow.c <- rep(NA,cols)
cihigh.c <- rep(NA, cols)
for (i in 1:cols) {
  mean.t[i] <- t.test(outcome.df[,i]~outcome.df$shut_off_planned)$estimate[2]
  cilow.t[i] <-  t.test(outcome.df[outcome.df$shut_off_planned==1,][i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df[outcome.df$shut_off_planned==1,][i])$conf.int[2]
  mean.c[i] <- t.test(outcome.df[,i]~outcome.df$shut_off_planned)$estimate[1]
  cilow.c[i] <-  t.test(outcome.df[outcome.df$shut_off_planned==0,][i])$conf.int[1]
  cihigh.c[i] <- t.test(outcome.df[outcome.df$shut_off_planned==0,][i])$conf.int[2]
  p.val[i] <- t.test(outcome.df[,i]~outcome.df$shut_off_planned)$p.value
}
means <- c(mean.t, mean.c)
cilow <- c(cilow.t,cilow.c)
cihigh <- c(cihigh.t, cihigh.c)
treat <- c(rep("Treat", cols), rep("Control", cols))
plotdata <- data.frame(rep(names(outcome.df)[1:ncol(outcome.df)-1],2) , means, cilow, cihigh, treat)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='already_backup', 1,
                              ifelse(outcome=='already_battery', 2,
                                     ifelse(outcome=='already_solar', 3,
                                            ifelse(outcome=='already_landscaping', 4,
                                                   ifelse(outcome=='already_ev', 6, NA))))))
plotdata$order <- cols+1-plotdata$order
plotdata$zone <- ifelse(plotdata$treat=="Treat", "Outage-exposed respondents", "Matched non-outage respondents")
plotdata <- arrange(plotdata, order)
plotdata$x <- sort(c((1:cols-.1),(1:cols+.1)))
outcomes2 <- c("Electric vehicle", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Backup generation")

prior.adaptation <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=zone))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=zone))+
  coord_flip() +
  theme_classic()+
  labs(x = "", y = "Proportion that already took adaptation behavior")+
  ylim(0,.7)+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  scale_color_manual(values=c('grey','black'))+
  theme(legend.position = "bottom", legend.title=element_blank())+
  theme(text = element_text(size=14)) #+
prior.adaptation
ggsave("figs/prior.adaptation.png", prior.adaptation, width=7, height=5)


###########################################################
##### TREATMENT DETERMINED BY LOCATION VS. SELFREPORT-- SI4 #####
############################################################

#select covariates 
X <- nosoc %>%
  dplyr::select(pid5, married, employed, non_english, kids_in_household, ideo, ed, female, smoke, age, ResponseId, sample, treat) %>%
  na.omit

#set.seed(12345)
#m.out.loc <- matchit(treat ~ pid5+ married+ employed+ non_english+ kids_in_household+ ideo+ ed+ female + smoke + age, data = X,
#                   pop.size=10000, method = "genetic", replace = T, estimand = "ATT") 

#mdtemp <- match.data(m.out.loc)
#export(mdtemp, "m.out.loc.xlsx")
mdtemp <- import("m.out.loc.xlsx")
m.data.loc <- left_join(subset(mdtemp, select = ResponseId), df, by="ResponseId") 


### adaptation analysis 

outcome.df <- m.data.loc %>% dplyr::select(battery, move_binary, solar_panels, ev, backup, extras, landscaping,treat) 
cols <- ncol(outcome.df)-1

mean.t = rep(NA, cols)
mean.c = rep(NA, cols)
p.val = rep(NA, cols)
cilow.t <- rep(NA, cols)
cihigh.t <- rep(NA, cols)
cilow.c <- rep(NA,cols)
cihigh.c <- rep(NA, cols)
for (i in 1:cols) {
  mean.t[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[2]
  cilow.t[i] <-  t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[2]
  mean.c[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[1]
  cilow.c[i] <-  t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[1]
  cihigh.c[i] <- t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[2]
  p.val[i] <- t.test(outcome.df[,i]~outcome.df$treat)$p.value
}
means <- c(mean.t, mean.c)
cilow <- c(cilow.t,cilow.c)
cihigh <- c(cihigh.t, cihigh.c)
treat <- c(rep("Treat", cols), rep("Control", cols))
plotdata <- data.frame(rep(names(outcome.df)[1:ncol(outcome.df)-1],2) , means, cilow, cihigh, treat)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='backup', 1,
                              ifelse(outcome=='extras', 2,
                                     ifelse(outcome=='battery', 3,
                                            ifelse(outcome=='solar_panels', 4, 
                                                   ifelse(outcome=='landscaping', 5,
                                                          ifelse(outcome=='move_binary', 6,
                                                                 ifelse(outcome=='ev', 7, NA))))))))
plotdata$order <- 8-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- sort(c((1:cols-.1),(1:cols+.1)))
plotdata$zone <- ifelse(plotdata$treat=="Treat", "Outage zone", "Matched non-outage zone")
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")

adaptation.plot.loc <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=zone))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=zone))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "Proportion intending to take adaptation behavior")+
  ylim(0,.7)+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  scale_color_manual(values=c('grey','black'))+
  theme(legend.title=element_blank())+
  ggplot2::annotate("text", x = 7.2, y = .18, label = "**")+
  ggplot2::annotate("text", x = 5.2, y = .10, label = "*")+
  ggplot2::annotate("text", x = 1.2, y = .22, label = "**")+
  theme(legend.position = "bottom", legend.title=element_blank())
adaptation.plot.loc
ggsave("figs/adaptation.plot.loc.png", adaptation.plot.loc, width=7, height=5)

#p-vals
options(scipen=100)
t.test(m.data.loc$backup ~m.data.loc$treat)$p.value
t.test(m.data.loc$battery ~m.data.loc$treat)$p.value
t.test(m.data.loc$ev ~m.data.loc$treat)$p.value


### utility trust
m.data.loc$utility_trust_binary2 <- 1 - m.data.loc$utility_trust_binary
outcome.df <- m.data.loc %>% dplyr::select(utility_trust_binary2, utility_responsible_binary, pge_liable_binary, pge_role_binary, treat) #utility responsibility not as good b/c makes less sense for SoCal 
cols <- ncol(outcome.df)-1

mean.t = rep(NA, cols)
mean.c = rep(NA, cols)
p.val = rep(NA, cols)
cilow.t <- rep(NA, cols)
cihigh.t <- rep(NA, cols)
cilow.c <- rep(NA,cols)
cihigh.c <- rep(NA, cols)
for (i in 1:cols) {
  mean.t[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[2]
  cilow.t[i] <-  t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[1]
  cihigh.t[i] <- t.test(outcome.df[outcome.df$treat==1,][i])$conf.int[2]
  mean.c[i] <- t.test(outcome.df[,i]~outcome.df$treat)$estimate[1]
  cilow.c[i] <-  t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[1]
  cihigh.c[i] <- t.test(outcome.df[outcome.df$treat==0,][i])$conf.int[2]
  p.val[i] <- t.test(outcome.df[,i]~outcome.df$treat)$p.value
}
means <- c(mean.t, mean.c)
cilow <- c(cilow.t,cilow.c)
cihigh <- c(cihigh.t, cihigh.c)
treat <- c(rep("Treat", cols), rep("Control", cols))
plotdata <- data.frame(rep(names(outcome.df)[1:ncol(outcome.df)-1],2) , means, cilow, cihigh, treat)
names(plotdata)[1] <- "outcome"

plotdata$order <- with(plotdata,
                       ifelse(outcome=='utility_trust_binary2', 1,
                              ifelse(outcome=='utility_responsible_binary', 2,
                                     ifelse(outcome=='pge_liable_binary', 3, 
                                            ifelse(outcome=='pge_role_binary', 4, NA)))))

plotdata$order <- (cols+1)-plotdata$order
plotdata <- arrange(plotdata, order)
plotdata$x <- c(.9,1.1,1.9,2.1,2.9,3.1, 3.9, 4.1)
plotdata$zone <- ifelse(plotdata$treat=="Treat", "Outage zone", "Matched non-outage zone")
outcomes2 <- c( "PG&E major restructuring", "PG&E liable for wildfires (strongly agree)" , "Utility completely responsible for outages",  "Completely distrust their utility")

utilities.loc <- ggplot(plotdata, aes(x = x))+
  geom_point(size=3, aes(y=means, color=zone))+
  geom_linerange(aes(ymin=cilow,ymax=cihigh, color=zone))+
  coord_flip()+
  theme_classic()+
  labs(x = "", y = "Proportion taking position")+
  scale_x_continuous(breaks = 1:cols, labels = outcomes2)+ 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))+
  scale_color_manual(values=c('grey','black'))+
  theme(legend.title=element_blank())+
  ggplot2::annotate("text", x = 4.2, y = .4, label = "**")+
  ggplot2::annotate("text", x = 3.2, y = .65, label = "**")+
  theme(text = element_text(size=14)) +
  theme(legend.position = "bottom", legend.title=element_blank())
utilities.loc
ggsave("figs/utilities.loc.png", utilities.loc, width=7, height=4)

#p-vals
options(scipen=100)
t.test(m.data.loc$utility_trust_binary ~m.data.loc$treat)$p.value
t.test(m.data.loc$utility_responsible ~m.data.loc$treat)$p.value



### 2sls adaptation 
battery <- ivreg(battery ~ shut_off_planned | treat , data = m.data.loc)
backup <- ivreg(backup ~ shut_off_planned | treat , data = m.data.loc)
move_binary <- ivreg(move_binary ~ shut_off_planned | treat , data = m.data.loc)
solar_panels <- ivreg(solar_panels ~ shut_off_planned | treat , data = m.data.loc)
ev <- ivreg(ev ~ shut_off_planned | treat , data = m.data.loc)
extras <- ivreg(extras ~ shut_off_planned | treat , data = m.data.loc)
landscaping <- ivreg(landscaping ~ shut_off_planned | treat , data = m.data.loc)

coef <- c(ev$coefficients[2], move_binary$coefficients[2], landscaping$coefficients[2], solar_panels$coefficients[2],
          battery$coefficients[2], extras$coefficients[2], backup$coefficients[2])
se <- c(robust.se(ev)[2,2], robust.se(move_binary)[2,2], robust.se(landscaping)[2,2], robust.se(solar_panels)[2,2], 
        robust.se(battery)[2,2], robust.se(extras)[2,2], robust.se(backup)[2,2])
order <- c(1:7)
pdata <- data.frame(coef, se, order)
outcomes2 <- c("Electric vehicle", "Move", "Wildfire-robust landscaping", "Solar panels", "Home battery", "Extra food/water", "Backup generation")
cols <- length(outcomes2)

iv.adapt <- ggplot(pdata, aes(x = order, y = coef))+
  geom_point(size=3)+
  geom_linerange(aes(ymin=coef-1.96*se, ymax=coef+1.96*se), linetype = 1, size = 1)+
  geom_linerange(aes(ymin=coef-1.645*se, ymax=coef+1.645*se), linetype = 1, size = 1.5)+  coord_flip()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme_classic()+
  labs(x = "", y = "Estimated effect of outage on adaptation intentions")+
  ylim(-.5,.5)+
  scale_x_continuous(breaks = 1:cols, labels=outcomes2)+ #, labels = outcomes2
  theme(text = element_text(size=14))+
  scale_color_manual(values=c('black','grey'))+
  theme(legend.position = "bottom", legend.title=element_blank())
iv.adapt
ggsave("figs/iv.adapt.png", iv.adapt, width=7, height=5)

### utility trust
m.data.loc$utility_trust_binary2 <- 1 - m.data.loc$utility_trust_binary

trust <- ivreg(utility_trust_binary2 ~ shut_off_planned | treat, data = m.data.loc)
responsible <- ivreg(utility_responsible_binary ~ shut_off_planned | treat, data = m.data.loc)
liable <- ivreg(pge_liable_binary ~ shut_off_planned | treat, data = m.data.loc)
role <- ivreg(pge_role_binary ~ shut_off_planned | treat, data = m.data.loc)

coef <- c(role$coefficients[2], liable$coefficients[2], responsible$coefficients[2], trust$coefficients[2])
se <- c(robust.se(role)[2,2], robust.se(liable)[2,2], robust.se(responsible)[2,2], robust.se(trust)[2,2])
order <- c(1:4)
pdata <- data.frame(coef, se, order)
outcomes2 <- c( "PG&E major restructuring", "PG&E liable for wildfires (strongly agree)" , "Utility completely responsible for outages",  "Completely distrust their utility")
cols <- length(outcomes2)


iv_utility <- ggplot(pdata, aes(x = order, y = coef))+
  geom_point(size=3)+
  geom_linerange(aes(ymin=coef-1.96*se, ymax=coef+1.96*se), linetype = 1, size = 1)+
  geom_linerange(aes(ymin=coef-1.645*se, ymax=coef+1.645*se), linetype = 1, size = 1.5)+  coord_flip()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme_classic()+
  labs(x = "", y = "Estimated effect of outage on PG&E attitudes")+
  ylim(-.5,.5)+
  scale_x_continuous(breaks = 1:cols, labels=outcomes2)+ #, labels = outcomes2
  theme(text = element_text(size=14))+
  scale_color_manual(values=c('black','grey'))+
  theme(legend.position = "bottom", legend.title=element_blank())
iv_utility
ggsave("figs/iv.utility.png", iv_utility, width=7, height=5)





######################################################
##### WILLINGNESS TO PAY ANALYSIS ##########
#################################################

#new sample coding
df$sample2 <- as.factor(with(df,
                   ifelse(sample %in% c("inner.all","inner.0.1"), "inner",
                          ifelse(sample=="socal", "socal", "outer"))))

#days w/o electricity
days.without <- sbchoice(no_electricity_reduce ~ sample2 | no_electricity_condition, dist = "logistic", data = subset(df, sample!="socal"))
summary(days.without)

#monthly surcharge to avoid future outages
surcharge <- sbchoice(surcharge_reduce ~ sample2 | surcharge_condition, dist = "logistic", data = subset(df, sample!="socal"))
summary(surcharge)

#monthly fee to bury
bury <- sbchoice(bury_reduce ~ sample2 | bury_condition, dist = "logistic", data = subset(df, sample!="socal"))
summary(bury)

